#!/usr/bin/python2.5
"""

parse_ncbi_exon_data.pl
    --exon_input_file PATH
    --cdna_input_file PATH
    --prot_input_file PATH

"""

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   options        


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

from optparse import OptionParser
import sys, os
import os.path
import StringIO








# add self to search path for testing
if __name__ == '__main__':
    exe_path = os.path.split(os.path.abspath(sys.argv[0]))[0]
    sys.path.append(os.path.abspath(os.path.join(exe_path,"..", "python_modules")))
    myname = os.path.split(sys.argv[0])[1]
    myname = os.path.splitext(myname)[0];
else:
    myname = __name__

import general_util
from general_util import die_error,prepare_directory_for_file
from custom_log import custom_log
from t_ncbi_gene_coding_structure import t_gene, t_locus, t_transcript
import t_ncbi_gene_coding_structure
from get_genomic_sequence import get_genomic_sequence_coordinates_from_index_file,get_genomic_sequence




parser = OptionParser(version="%prog 1.0")


#
#   exons
# 

parser.add_option("-l", "--group_label", dest = "group_labels",
                    action="append", 
                    metavar="LABEL", 
                    type="string",
                    help="The protocol of interest, e.g. reference, build36 etc.")
parser.add_option("-I", "--input_directory", dest="input_directory",
                  metavar="PATH", 
                  type="string",
                  help="path of input directory (assumes cdna.orig.fa and prot.orig.fa for "
                    "input files)")
parser.add_option("-O", "--output_directory", dest="output_directory",
                  metavar="PATH", 
                  type="string",
                  help="path of output directory (assumes cdna.fa and prot.fa, "+
                        "and exon_loci.data for the output files)")
parser.add_option("-W", "--working_directory", dest="working_directory",
                  metavar="PATH", 
                  type="string",
                  help="path of working directory (assumes exonic_sequence.cache, "
                       "exonic_loci.cache, parse.log)")          
parser.add_option("-N", "--name_prefix", dest="name_prefix",
                  metavar="PREFIX", 
                  type="string",
                  help="prefix for identifiers. Assumes genes, proteins and transcripts "+
                    "identifiers = PREFIX + G/P/T + '_'")          

# 
#   sequence input
# 
parser.add_option("-c", "--cdna_input_file", dest="cdna_input_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of cDNA file (rna.fa)")
parser.add_option("-p", "--prot_input_file", dest="prot_input_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of peptide file (protein.fa)")
parser.add_option("-e", "--exon_input_file", dest="exon_input_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of exon file (seq_gene.md)")

#   
#   sequence output
# 
parser.add_option("-C", "--cdna_output_file", dest="cdna_output_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of file which will hold the genomic CDS")
parser.add_option("-P", "--prot_output_file", dest="prot_output_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of file which will hold the translated genomic CDS")
parser.add_option("-E", "--exon_output_file", dest="exon_output_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of file to hold exonic loci [gene_id, "
                        "cdna_id, prot_id, chromosome, beg, end, strand] in 0-based")
parser.add_option("-R", "--refseq_output_file", dest="refseq_output_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of file to hold refseq associations[gene_id, "
                        "cdna_id, prot_id, cdna_refseq, prot_refseq] for sequences "
                        "where refseq and genomic sequence actually match!")
parser.add_option("-M", "--mapped_identifiers_output_file", dest="mapped_identifiers_output_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of file to hold gene_id, "
                        "cdna_id, prot_id")
parser.add_option("-Z", "--pseudos_with_cds_output_file", dest="pseudos_with_cds_output_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of file to hold list of pseudo genes with cds "
                        "[gene_id, reason]")

parser.add_option("-G", "--gis_for_cdna_prot_identifiers", dest="gis_for_cdna_prot_identifiers",
                    action="store_true", default=False,
                    help="Use GIs rather than refseq identifiers to construct cdna and prot identifiers.")


parser.add_option("--gene_accession_prefix", dest="gene_accession_prefix",
                  metavar="PREFIX", 
                  type="string",
                  help="Prefix will be prepended to GeneId to give gene identifier")
parser.add_option("--cdna_accession_prefix", dest="cdna_accession_prefix",
                  metavar="PREFIX", 
                  type="string",
                  help="Prefix will be prepended to cdna gi to give cdna sequence identifier")
parser.add_option("--prot_accession_prefix", dest="prot_accession_prefix",
                  metavar="PREFIX", 
                  type="string",
                  help="Prefix will be prepended to prot gi to give prot sequence identifier")

#
#   genomic sequence
# 
parser.add_option("-g", "--genomic_sequence_file", dest = "genomic_sequence_file",
                    metavar="FILE", 
                    type="string",
                    help="name and path of genomic sequence file (genomic.fa)")
parser.add_option("-i", "--genomic_index_file", dest = "genomic_index_file",
                    metavar="FILE", 
                    type="string",
                    help="name and path of genomic index file (genomic.index)")

#
#   caches
# 
parser.add_option("--exonic_sequence_cache_file", dest="exonic_sequence_cache_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of file used to cache the exonic sequence")
parser.add_option("--exonic_loci_cache", dest="exonic_loci_cache_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of file used to cache the exonic loci")


#
#   logger
# 
parser.add_option("-q", "--quiet", dest = "quiet",
                  action="store_true", default=False,
                  help="Do not echo to shell but only print to log.")

parser.add_option("-L", "--log_file", dest="log_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of log file")
parser.add_option("--log_parameters", dest="log_parameters",
                    action="store_true", default=False,
                    help="Print program parameters to log.")

#
#   debug dump
# 
parser.add_option("-1", "--mismatched_genes_dump_file", dest="mismatched_genes_dump_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of dump of genes where the genomic and "
                        "transcript are different")
parser.add_option("-2", "--disrupted_genes_dump_file", dest="stops_genes_dump_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of dump of genes with stops in the genomic cds")
parser.add_option("-3", "--all_genes_dump_file", dest="all_genes_dump_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of dump of all genes")
parser.add_option("-4", "--all_gene_types_dump_file", dest="all_gene_types_dump_file",
                  metavar="FILE", 
                  type="string",
                  help="name and path of dump of the types (CDS/pseudo etc.) of all genes "
                        "[gene_type, gene_id]")


parser.add_option("--debug", dest="debug_mode",
                action="store_true", default=False,
                help="Only iterates first 500000 lines and ignores cache.")
    

parameters = [  "exon_input_file",
                "group_labels",
                "cdna_input_file",
                "prot_input_file",
                "cdna_output_file",
                "prot_output_file",
                "exon_output_file",
                "refseq_output_file",
                "mapped_identifiers_output_file",
                "pseudos_with_cds_output_file",
                "prot_accession_prefix",
                "cdna_accession_prefix",
                "gene_accession_prefix",
                "genomic_sequence_file",
                "genomic_index_file",
                "exonic_loci_cache_file",
                "exonic_sequence_cache_file",
                "log_file",
                "stops_genes_dump_file"      ,
                "mismatched_genes_dump_file" ,
                "all_gene_types_dump_file",
                "all_genes_dump_file",
                "working_directory",
                "input_directory",
                "output_directory",
                "name_prefix",
                ]

# parameters default to none
for parameter in parameters:
    parser.set_defaults(**{parameter: None}) 

    
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#   DEBUG test parameters for ncbi homo
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
base_dir = "/net/cpp-mirror/databases/ensembl/homo_ncbi_36_3/"
parser.set_defaults(input_directory     = base_dir + "ncbi/")
parser.set_defaults(output_directory    = base_dir)
parser.set_defaults(working_directory   = base_dir + "parse_ncbi/")
parser.set_defaults(name_prefix         = "NCBIHSA")

parser.set_defaults(genomic_sequence_file       = base_dir + "dna/genomic.fa")
parser.set_defaults(genomic_index_file          = base_dir + "dna/genomic.index")
parser.set_defaults(all_gene_types_dump_file    = base_dir + "dump_all_genes.types")
#parser.set_defaults(group_labels = ["reference"])
parser.set_defaults(log_parameters = True)
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#   DEBUG test parameters for ncbi homo
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888





#   
#   get help string
#
f =StringIO.StringIO()
parser.print_help(f)
helpstr = f.getvalue()
(options, args) = parser.parse_args()

#options.group_labels = set(options.group_labels)

#_________________________________________________________________________________________
#
#   log file
# 
#       put here before helper options increase the number of defined options
#
#_________________________________________________________________________________________
if options.working_directory:
    if not options.log_file:
        options.log_file = os.path.join(options.working_directory, "parse.log")

# make sure there is a directory ready for output
prepare_directory_for_file(options.log_file)
status_log = custom_log(myname, options.log_file, options.quiet)
t_ncbi_gene_coding_structure.status_log = status_log
#status_log.no_new_line_on_exit()

#DEBUGG
if options.group_labels == None:
    options.group_labels = ["reference"]

#
#     Make a record of the program parameters so we can rerun things if necessary
#
from  get_option_strings import get_option_strings
if options.log_parameters:
    programme_option_parameters = myname + ".py \\\n" + \
                                " \\\n".join(get_option_strings (parser, options))
    status_log.info_unwrapped_to_log(programme_option_parameters + "\n")

#_________________________________________________________________________________________
#
#   Helper options
#
#_________________________________________________________________________________________
if options.output_directory:
    if not options.cdna_output_file:
        options.cdna_output_file = os.path.join(options.output_directory, "cdna.fa")
    if not options.prot_output_file:
        options.prot_output_file = os.path.join(options.output_directory, "prot.fa")
    if not options.exon_output_file:
        options.exon_output_file = os.path.join(options.output_directory, "exon_loci.data")
    if not options.refseq_output_file:
        options.refseq_output_file = os.path.join(options.output_directory, "refseq_xref.map")
    if not options.mapped_identifiers_output_file:
        options.mapped_identifiers_output_file = os.path.join(options.output_directory, "gene_prot_cdna_id.data")
    if not options.pseudos_with_cds_output_file:
        options.pseudos_with_cds_output_file = os.path.join(options.output_directory, "gene_ids.pseudos")
if options.input_directory:
    if not options.cdna_input_file:
        options.cdna_input_file = os.path.join(options.input_directory, "cdna.orig.fa")
    if not options.prot_input_file:
        options.prot_input_file = os.path.join(options.input_directory, "prot.orig.fa")
    if not options.exon_input_file:
        options.exon_input_file = os.path.join(options.input_directory, "seq_gene.md")

if options.working_directory:
    if not options.exonic_loci_cache_file:
        options.exonic_loci_cache_file = os.path.join(options.working_directory, 
                                                                "exonic_loci.cache")
    if not options.exonic_sequence_cache_file:
        options.exonic_sequence_cache_file = os.path.join(options.working_directory, 
                                                                "exonic_sequence.cache")
    if not options.log_file:
        options.log_file = os.path.join(options.working_directory, "parse.log")
if options.name_prefix:
    if not options.prot_accession_prefix:
        options.prot_accession_prefix = options.name_prefix + "P_"
    if not options.cdna_accession_prefix:
        options.cdna_accession_prefix = options.name_prefix + "T_"
    if not options.gene_accession_prefix:
        options.gene_accession_prefix = options.name_prefix + "G_"
        
if options.debug_mode:
    options.exonic_loci_cache_file = None
    options.exonic_sequence_cache_file = None
    
        

#   mandatory options
#  
mandatory_parameters =  [   "exon_input_file",
                            "group_labels",
                            "cdna_input_file",
                            "prot_input_file",
                            "cdna_output_file",
                            "prot_output_file",
                            "exon_output_file",
                            "prot_accession_prefix",
                            "cdna_accession_prefix",
                            "gene_accession_prefix",
                            "genomic_sequence_file",
                            "genomic_index_file",
                            ]


for parameter in mandatory_parameters:
    if options.__dict__[parameter] == None:
        die_error("Please specify a file in --%s.\n\n" % parameter + helpstr)






#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   imports        


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

import glob
import shutil
import cPickle
import re
from collections import defaultdict
from general_util import die_error
from sequence_locus import *









#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   globals


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888










#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Functions        


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#_________________________________________________________________________________________
# 
#   print dots to show progress
#
#_________________________________________________________________________________________
def progress_dots (cnt_genes_dots):
    if cnt_genes_dots[0] % cnt_genes_dots[2] == 0:
        cnt_genes_dots[1] += 1
        if cnt_genes_dots[1] % 50 ==0:
            sys.stderr.write("\n    ")
        sys.stderr.write(".")

count_PAR_gene_duplicates = 0
count_buggy_gene_duplicates = 0
buggy_gene_duplicates = list()
PAR_gene_duplicates = list()
#_________________________________________________________________________________________
# 
#   get_exonic_structure_from_cache_file
#
#_________________________________________________________________________________________
def get_exonic_structure_from_cache_file ():
    """
        if exonic structure already exists in cache, use that
        else 
            parse seq_gene.md  

    """
    global count_PAR_gene_duplicates
    global count_buggy_gene_duplicates
    global buggy_gene_duplicates
    global PAR_gene_duplicates

    #   get genomic sequence per coding exon
    # 
    if options.exonic_loci_cache_file and os.path.exists(options.exonic_loci_cache_file):

        status_log.info("Loading exonic structure from \n[%s]" %
                            options.exonic_loci_cache_file)
        cache_file = open(options.exonic_loci_cache_file, 'rb')
        exonic_structures = cPickle.load(cache_file)
        gene_id_versions = cPickle.load(cache_file)
        count_PAR_gene_duplicates = cPickle.load(cache_file)
        count_buggy_gene_duplicates = cPickle.load(cache_file)
        buggy_gene_duplicates = cPickle.load(cache_file)
        PAR_gene_duplicates = cPickle.load(cache_file)
        cache_file.close()
    else:

        #   parse seq_gene.md  
        # 
        #           => exonic_structures[group_label][gene_id] = [t_gene, t_gene, t_gene]
        #
        exon_file = open(options.exon_input_file)
        ignore_headers_on_first_line = exon_file.readline()
        exonic_structures = defaultdict(dict)
        #
        #   holds gene_id versions identifying separate gene records with the same gene_id
        gene_id_versions = defaultdict(lambda : defaultdict(list))
        parse_exon_structure_file (exonic_structures, gene_id_versions, exon_file, options.group_labels)


        # 
        #   don't use defaultdict because can't pickle 
        #
        picklable_gene_id_versions = dict()
        for group_label in gene_id_versions:
            picklable_gene_id_versions[group_label] = dict()
            for original_id, versioned_gene_ids in gene_id_versions[group_label].iteritems():
                picklable_gene_id_versions[group_label][original_id] = versioned_gene_ids
        gene_id_versions = picklable_gene_id_versions

        #
        #   save to cache 
        #
        if options.exonic_loci_cache_file:
            status_log.info("Caching exonic structure to \n[%s]" %
                                options.exonic_loci_cache_file)
            # make sure there is a directory ready for output
            prepare_directory_for_file(options.exonic_loci_cache_file)
            cache_file = open(options.exonic_loci_cache_file, 'wb')
            cPickle.dump(exonic_structures, cache_file, 1)
            cPickle.dump(gene_id_versions, cache_file, 1)
            cPickle.dump(count_PAR_gene_duplicates  , cache_file, 1)
            cPickle.dump(count_buggy_gene_duplicates, cache_file, 1)
            cPickle.dump(buggy_gene_duplicates      , cache_file, 1)
            cPickle.dump(PAR_gene_duplicates      , cache_file, 1)
            cache_file.close()

    return exonic_structures, gene_id_versions



#_________________________________________________________________________________________
# 
#   get_exonic_sequence_from_cache_file
#
#_________________________________________________________________________________________
def get_exonic_sequence_from_cache_file (genes):
    """
        if exonic sequence already exists in cache, use that
        
    """
    #   get genomic sequence per coding exon
    # 
    if options.exonic_sequence_cache_file and os.path.exists(options.exonic_sequence_cache_file):

        status_log.info("Loading exonic sequence from \n[%s]" %
                            options.exonic_sequence_cache_file)
        cache_file = open(options.exonic_sequence_cache_file, 'rb')
        exonic_coding_sequence = cPickle.load(cache_file)
        coding_exonic_sequence_errors = cPickle.load(cache_file)
        cache_file.close()
    else:
        index_file = open(options.genomic_index_file)               
        seq_coordinates = get_genomic_sequence_coordinates_from_index_file (index_file)
        exonic_coding_loci =  get_exonic_coding_loci (genes)
        status_log.info("Loading exonic sequence for %d genes %d exons" % 
                                        (len(genes), len(exonic_coding_loci)))
        (exonic_coding_sequence,
        coding_exonic_sequence_errors) = get_exonic_coding_sequence (exonic_coding_loci, 
                                                            seq_coordinates, 
                                                            options.genomic_sequence_file)
        #
        #   DEBUGG sanity check 
        #
        for locus in exonic_coding_loci:
            if not exonic_coding_sequence.has_key(locus):
                if not coding_exonic_sequence_errors.has_key(locus):
                    sys.stderr.write(str(locus))
                    raise RuntimeError, locus

        #
        #   save to cache 
        #
        if options.exonic_sequence_cache_file:
            status_log.info("Caching exonic sequence to \n[%s]" %
                                options.exonic_sequence_cache_file)
            # make sure there is a directory ready for output
            prepare_directory_for_file(options.exonic_sequence_cache_file)
            cache_file = open(options.exonic_sequence_cache_file, 'wb')
            cPickle.dump(exonic_coding_sequence, cache_file, 1)
            cPickle.dump(coding_exonic_sequence_errors, cache_file, 1)
            cache_file.close()

    return exonic_coding_sequence, coding_exonic_sequence_errors
#_________________________________________________________________________________________
# 
#   get_exonic_coding_loci
#
#_________________________________________________________________________________________
def get_exonic_coding_loci (genes):
    status_log.debug("Get loci for each coding exon...")
    coding_loci = set()
    for g in genes:
        coding_loci.update(g.get_coding_loci_before_finalise_data())
    coding_loci = list(coding_loci)
    coding_loci.sort()
    return coding_loci


#_________________________________________________________________________________________
# 
#   get_exonic_coding_sequence
#
#_________________________________________________________________________________________
import copy
def get_exonic_coding_sequence (coding_loci, seq_coordinates, seq_file):

    status_log.debug("Get sequence for each coding exon...")

    coding_exonic_sequence = dict()
    coding_exonic_sequence_errors = dict()

    # progress dots begin
    if not options.quiet:
        sys.stderr.write("    [")
    cnt_genes_dots = [0,0, 1]
    cnt_genes_dots[2] = int(len(coding_loci) / 100) + 1

    for locus in coding_loci:
        try:
            exon_seq, \
            chromosome = get_genomic_sequence (locus.chromosome,
                                                    locus.beg, locus.len,
                                                    seq_coordinates, seq_file)                
            coding_exonic_sequence[locus] = [exon_seq, chromosome]
            # use new chromosome as well
            if chromosome != None:
                new_locus = copy.copy(locus)
                new_locus.chromosome = chromosome
                coding_exonic_sequence[new_locus] = [exon_seq, None]
    
            # progress dots
            cnt_genes_dots[0] += 1
            if not options.quiet and (cnt_genes_dots[0] % cnt_genes_dots[2]) == 0:
                progress_dots(cnt_genes_dots)
                
        except RuntimeError, e:
            status_log.warning (e.message)
            coding_exonic_sequence_errors[locus] = e.message

    # progress dots end
    if not options.quiet:
        sys.stderr.write("]\n")


    return coding_exonic_sequence, coding_exonic_sequence_errors

#_________________________________________________________________________________________
# 
#   get_group_labels_genes
#
#_________________________________________________________________________________________
def get_group_labels_genes (exonic_structures, group_labels):
    genes = list()
    for group_label in group_labels:
        for gene_id in exonic_structures[group_label].iterkeys():
            gene = exonic_structures[group_label][gene_id]
            genes.append(gene)
    return genes


#_________________________________________________________________________________________
# 
#   add_genomic_coding_sequence
#
#_________________________________________________________________________________________
def add_genomic_coding_sequence (genes, exonic_coding_sequence, coding_exonic_sequence_errors):
    """
        go through all genes and add cds genomic sequence
    """
    status_log.debug("Add genomic sequence...")

    # progress dots begin
    if not options.quiet:
        sys.stderr.write("    [")
    cnt_genes_dots = [0,0, 200]

    for gene in genes:
        gene.add_genomic_coding_sequence(exonic_coding_sequence, coding_exonic_sequence_errors)

        # progress dots
        cnt_genes_dots[0] += 1
        if not options.quiet and (cnt_genes_dots[0]) % cnt_genes_dots[2] == 0:
            progress_dots(cnt_genes_dots)

    # progress dots end
    if not options.quiet:
        sys.stderr.write("]\n")
            
            


#_________________________________________________________________________________________
# 
#   get_coding_sequence_statistics
#
#_________________________________________________________________________________________
def print_stats (stats, evidence):

    categories = sorted(stats["all"].keys())
    chromosomes = sorted(stats["chromosome"].keys(), cmp = lambda x,y: cmp(len(x),len(y)) or cmp(x,y))
    # 
    # DEBUG
    #   only print out summaries for chromosome 1, 2, x, y
    chromosomes = [c for c in chromosomes if c in ["1", "2", "x", "X", "y", "Y", "MT"]]
    
    
    strands = ('+', '-')
    
    all_cat = " 1 genes with cds"
    # all genes    
    status_log.info("All genes")
    for cat in categories:
        if stats["all"][cat]:
            status_log.info("  %-40s = %8d (%5.1f%%)" % (cat, stats["all"][cat], 
                                                        stats["all"][cat] * 100.0 /
                                                        stats["all"][all_cat]))
    # by strand    
    status_log.info("By strand")
    for strand in (True, False):
        for cat in categories:
            if stats["strand"][strand][cat]:
                status_log.info("  %-40s[%s] = %8d (%5.1f%%)" % (cat, strands[strand], 
                                                        stats["strand"][strand][cat], 
                                                        stats["strand"][strand][cat] * 100.0 /
                                                        stats["strand"][strand][all_cat]))
    # by chromosome    
    status_log.info("By chromosome")
    for chrm in chromosomes:
        if not chrm in stats["chromosome"]:
            continue
            
        for cat in categories:
            if stats["chromosome"][chrm][cat]:
                status_log.info("  %-40s[%s] = %8d (%5.1f%%)" % (cat, chrm, 
                                                        stats["chromosome"][chrm][cat], 
                                                        stats["chromosome"][chrm][cat] * 100.0 /
                                                        stats["chromosome"][chrm][all_cat]))
            
    for cat in ["matching", "mismatched", "stops"]:
        status_log.info("%s evidence:" % cat)
        for e in sorted(evidence[cat].keys()):
            status_log.info("  %-30s = %-10d" %(e, evidence[cat][e]))
    
def inc_stats (stats, strand, chrm, category):
    stats["strand"][strand][category] += 1
    stats["chromosome"][chrm][category] += 1
    stats["all"][category] += 1

def get_coding_sequence_statistics (genes):
    
    """
        go through all genes and summarise stats
    """
    
    
    
    status_log.debug("get_coding_sequence_statistics...")

    stats = dict()
    stats["strand"]     = defaultdict(lambda: defaultdict(lambda:0))
    stats["chromosome"] = defaultdict(lambda: defaultdict(lambda:0))
    stats["all"]        = defaultdict(lambda:0)
    
    evidence = dict()
    for t in ["matching", "mismatched", "stops", "tolerated_stops"]:
        evidence[t] = defaultdict(lambda:0)

    for gene in genes:
        
        #
        #   only include genes with cds 
        #
        if not gene.cnt_cds_transcripts:
            continue
            
            
        strand = gene.gene_locus.strand
        chrm = gene.gene_locus.chromosome
        inc_stats(stats, strand, chrm, " 1 genes with cds")

        # coding gene
        if gene.gene_type == "coding":
            inc_stats(stats, strand, chrm, " 2 Coding")
        else:
            inc_stats(stats, strand, chrm, " 3 " + gene.gene_type)
            
        if gene.cnt_known_genomic_stops:
            inc_stats(stats, strand, chrm, " 4 Transcripts with known stops")

        if gene.cnt_missing_prot_refseq:
            inc_stats(stats, strand, chrm, " 5 Transcripts with missing prot refseq")

        if gene.cnt_mismatched_genomic:
            inc_stats(stats, strand, chrm, " 6 Transcripts not matching genomic")
            
        if gene.cnt_genomic_stops:
            inc_stats(stats, strand, chrm, " 7 Transcripts with stops")
    
        if gene.missing_genomic_sequence:
            inc_stats(stats, strand, chrm, " 8 Missing genomic sequence")
            
        if gene.all_mismatched_genomic:
            inc_stats(stats, strand, chrm, " 9 all refseq != genomic")
        if gene.all_genomic_stops:
            inc_stats(stats, strand, chrm, "10 all genomic stops")
        if gene.all_genomic_stops and gene.all_mismatched_genomic:
            inc_stats(stats, strand, chrm, "11 all stops and mismatches")
        
        #
        #   collect evidence used for different categories
        #     
        for evidence_type, categories in gene.transcript_evidence.iteritems():
            for cat in categories:
                evidence[evidence_type][cat] += 1


    #
    #   statistics gathered in strands
    #       but code seems to be ok: ignore strand 
    #
    print_stats (stats, evidence)
    
#_________________________________________________________________________________________
# 
#   get_genes_with_stops_mismatches
#
#_________________________________________________________________________________________
def get_genes_with_stops_mismatches (genes):

    """
        go through all genes and get genes with stops and mismatches
    """

    stops_genes         = list()
    mismatched_genes    = list()


    status_log.debug("get_genes_with_stops_mismatches...")

    for gene in genes:
        if not gene.cnt_cds_transcripts:
            continue
        if gene.all_mismatched_genomic:
            mismatched_genes.append(gene)
        if gene.all_genomic_stops:
            stops_genes.append(gene)

    status_log.debug("%d genes with stops..." % len(stops_genes))
    status_log.debug("%d genes with mismatches..." % len(mismatched_genes))
    return stops_genes, mismatched_genes

#_________________________________________________________________________________________
# 
#   output_duplicate_coding_genes
#
#_________________________________________________________________________________________
def output_duplicate_coding_genes (gene_id_versions, exonic_structures):
    status_log.debug("Output duplicate coding genes...")
    XYgenes = list()
    non_XYgenes = list()
    cnt_non_XYgene_duplications = defaultdict(lambda: 0)
    
    
    for group_label in gene_id_versions:
        for original_gene_id in gene_id_versions[group_label].keys():
            PARxygene = False
            versioned_gene_ids = gene_id_versions[group_label][original_gene_id]
            cnt_coding_genes = 0
            for gene_id in versioned_gene_ids:
                gene = exonic_structures[group_label][gene_id]
                if gene.gene_type in set( ["coding", "coding sans ids"]):
                    cnt_coding_genes += 1
                    if gene.gene_id[-1] in "XYxy":
                        XYgenes.append("    %-10s : %-10s : %-10s : %-10s" % (gene.gene_type, 
                                                        gene.gene_id, 
                                                        group_label,
                                                        gene.gene_name))
                        PARxygene = True
                    else:
                        non_XYgenes.append(gene)
            if not PARxygene and cnt_coding_genes:
                cnt_non_XYgene_duplications[cnt_coding_genes] += 1

    status_log.info("%d PAR x/y chromosome genes" % (len(XYgenes) / 2))
    for g in XYgenes:
        status_log.info(g)
        
    for dupl_count in sorted(cnt_non_XYgene_duplications.keys()):
        status_log.info("%d other genes with %d copies" % (cnt_non_XYgene_duplications[dupl_count],
                                                            dupl_count))
    for g in non_XYgenes:
        status_log.info(str(g))


#_________________________________________________________________________________________
# 
#   output gene counts by type
#
#_________________________________________________________________________________________
def output_gene_counts_by_type (exonic_structures):
    """
 
    count "rna", "pseudo", "pseudo sans cds", "coding sans ids", "coding"
       for each group label
  

    """
    status_log.debug("Output gene counts by type...")
    gene_type_cnts = defaultdict(lambda: defaultdict(lambda:0))
    for group_label in exonic_structures:
        for gene in exonic_structures[group_label].itervalues():
            gene_type_cnts[group_label][gene.gene_type] += 1

    for group_label in gene_type_cnts:
        status_log.info(group_label)
        for gene_type in gene_type_cnts[group_label]:
            status_log.info("    %-20s = %d" % (gene_type, gene_type_cnts[group_label][gene_type]))


            
            
            
            
            
#_________________________________________________________________________________________
# 
#   get_refseq_ids_to_transcript
#
#_________________________________________________________________________________________
def get_refseq_ids_to_transcript (genes):
    """
        go through all transcripts in all genes and save the relationships between
            cdna_refseq_id   => transcript data structure
            prot_refseq_id   => transcript data structure
    """
    status_log.debug("Look up transcript and peptide identifiers...")

    prot_refseq_id_to_transcript =defaultdict(list)
    cdna_refseq_id_to_transcript =defaultdict(list)
        
    for gene in genes:
        for cdna_refseq_id, transcript in gene.transcripts.iteritems():
            prot_refseq_id = transcript.prot_refseq_id
            if prot_refseq_id != None and prot_refseq_id != '-':
                prot_refseq_id_to_transcript[prot_refseq_id].append(transcript)
            if cdna_refseq_id != None and cdna_refseq_id != '-':
                cdna_refseq_id_to_transcript[cdna_refseq_id].append(transcript)
            
            
    return prot_refseq_id_to_transcript, cdna_refseq_id_to_transcript
            
            
#_________________________________________________________________________________________
# 
#   get_gis_to_transcript
#
#_________________________________________________________________________________________
def get_gis_to_transcript (genes):
    """
        go through all transcripts in all genes and save the relationships between
            cdna_gi   => transcript data structure
            prot_gi   => transcript data structure
    """
    status_log.debug("Look up transcript and peptide GIs...")

    prot_gi_to_transcript =defaultdict(list)
    cdna_gi_to_transcript =defaultdict(list)

    for gene in genes:
        for transcript in gene.transcripts.itervalues():
            prot_gi = transcript.prot_gi
            cdna_gi = transcript.cdna_gi
            if prot_gi != None and prot_gi != '-':
                prot_gi_to_transcript[prot_gi].append(transcript)
            if cdna_gi != None and cdna_gi != '-':
                cdna_gi_to_transcript[cdna_gi].append(transcript)


    return prot_gi_to_transcript, cdna_gi_to_transcript


            
#_________________________________________________________________________________________
# 
#   resolve_duplicate_gis
#
#_________________________________________________________________________________________
def resolve_duplicate_gis (gi_to_transcript, gi_field_name):
    """
        go through all transcripts with duplicate gis and add version numbers
    """
    for gi, transcripts in gi_to_transcript.iteritems():
        if len(transcripts) == 1:
            continue
        
        for i, t in enumerate(transcripts):
            t.__dict__[gi_field_name] = "%sv%d" % (gi, i)



            
            




#_________________________________________________________________________________________
# 
#   dump_and_die
#
#       print out offending line parsed into separate fields and exit programme
#_________________________________________________________________________________________
def dump_and_die (reason, line_num, line):
    field_strs = ["tax_id", "chrm", "chr_start", "chr_stop", "strand", "contig", 
              "ctg_start", "ctg_stop", "ctg_orient", "feature_name", "gene_id", 
              "feature_type", "group_label", "transcript", "evidence_code"]
    fields = line.strip().split('\t')
    results = list()
    for i, field_val in enumerate(fields):
        results.append("    %-15s = [%s]\n" % (field_strs[i], field_val))
    status_log.die_error(reason + " on line %d\n%s" % 
                            (line_num, "".join(results)))


    
    

#_________________________________________________________________________________________
# 
#   new_gene_progress_dots
# 
#       helper function
#
#_________________________________________________________________________________________
def new_gene_progress_dots (cnt_genes_dots):
    """
    print progress dots for new genes
    """
    # progress dots
    cnt_genes_dots[0] += 1
    if not options.quiet and (cnt_genes_dots[0]) % cnt_genes_dots[2] == 0:
        progress_dots(cnt_genes_dots)


#_________________________________________________________________________________________
# 
#   create_duplicate_gene_record
#
#_________________________________________________________________________________________
def create_duplicate_gene_record (exonic_structures, gene_id_versions, 
                                    group_label, gene_id, locus,
                                    cnt_genes_dots):
    """
    
        When separate genes have the same gene_id, give each a separate gene_id with
            a different suffix (.xy for PAR, .1234 otherwise)
    
    """
    global count_buggy_gene_duplicates
    global count_PAR_gene_duplicates
    global buggy_gene_duplicates
    global PAR_gene_duplicates
    original_gene = exonic_structures[group_label][gene_id]
    
    new_gene_id_suffix = ""

    # 
    # second version of this gene
    #
    if not gene_id in gene_id_versions[group_label]:

        #
        #   genes on X and Y of PAR
        if (original_gene.gene_locus != None                and 
            original_gene.gene_locus.chromosome in "XYxy"   and
            locus.chromosome in "XYxy"                      and
            locus.chromosome != original_gene.gene_locus.chromosome):

            # save each version with a gene_id of its chromosome as subscript
            original_gene.gene_id_suffix  = ":" + original_gene.gene_locus.chromosome
            original_gene.gene_id        += original_gene.gene_id_suffix
            new_gene_id_suffix = ":" + locus.chromosome
            new_gene_id = gene_id + new_gene_id_suffix
            count_PAR_gene_duplicates += 1
            PAR_gene_duplicates.append(gene_id)
        else:
            # save each with version number as suffix
            original_gene.gene_id_suffix  = ":1"
            original_gene.gene_id        += ":1"
            new_gene_id_suffix = ":2"
            new_gene_id = gene_id + ":2"
            count_buggy_gene_duplicates += 1
            buggy_gene_duplicates.append(gene_id)

        #
        #   move original gene to correct place based on its new gene_id
        exonic_structures[group_label][original_gene.gene_id] = original_gene
        del exonic_structures[group_label][gene_id]
        gene_id_versions[group_label][gene_id].append(original_gene.gene_id)

    else:
        #
        #   nth version of this gene where n > 2
        #
        version_num = len (gene_id_versions[group_label][gene_id]) + 1
        new_gene_id_suffix = ":" + str(version_num)
        new_gene_id = gene_id + new_gene_id_suffix

    gene_id_versions[group_label][gene_id].append(new_gene_id)

    new_gene = t_gene(new_gene_id, group_label, locus, new_gene_id_suffix)
    exonic_structures[group_label][new_gene_id] = new_gene


    # progress dots
    new_gene_progress_dots (cnt_genes_dots)

    return new_gene
        
        
        
        
#_________________________________________________________________________________________
# 
#   create_or_retrieve_gene
#
#_________________________________________________________________________________________
def create_or_retrieve_gene (feature_type,
                                exonic_structures, gene_id_versions, 
                                group_label, gene_id, locus, cnt_genes_dots):
    if feature_type != "GENE":

        #
        #   1) Has duplicate (Multiple gene_ids map to same gene)
        # 
        #       return most recently encountered gene record
        #  
        if gene_id_versions[group_label].has_key(gene_id):
            gene_id = gene_id_versions[group_label][gene_id][-1]
            prexisting_gene =  exonic_structures[group_label][gene_id]


        #
        #   2) newly encountered gene_id: 
        #           create new gene record
        #  
        elif not gene_id in exonic_structures[group_label]:
            exonic_structures[group_label][gene_id] = t_gene(gene_id, group_label, locus)

            # progress dots
            new_gene_progress_dots (cnt_genes_dots)


            return exonic_structures[group_label][gene_id]
            
        #
        #   3) Existing gene
        #  
        else:
            prexisting_gene = exonic_structures[group_label][gene_id]

            
        #
        #   create new gene if on different chromosome / strand as existing gene
        #       even if not new gene record 
        #
        if (prexisting_gene.gene_locus.chromosome != locus.chromosome or 
            prexisting_gene.gene_locus.strand != locus.strand):
            return create_duplicate_gene_record (exonic_structures, gene_id_versions, 
                                                group_label, gene_id, locus, cnt_genes_dots)
        else:
            return prexisting_gene

        
        
    #------------------------------------------------------------------------------------
    # 
    # feature_type == "GENE"
    #
    if (gene_id in exonic_structures[group_label]):

        #
        #   oops different gene structure entry with same gene_id
        #
        if ( exonic_structures[group_label][gene_id].gene_name != None):

            return create_duplicate_gene_record (exonic_structures, gene_id_versions, 
                                                    group_label, gene_id, locus,
                                                    cnt_genes_dots)


        #
        #   otherwise just use this gene
        else:
            return exonic_structures[group_label][gene_id]
        
            
        
    #
    #   never encountered this gene_id before 
    #       make new gene structure
    #
    else:            
        # progress dots
        new_gene_progress_dots (cnt_genes_dots)

        exonic_structures[group_label][gene_id] = t_gene(gene_id, group_label, locus)
        return exonic_structures[group_label][gene_id]
        
    

#_________________________________________________________________________________________
# 
#   set_gene_and_transcript_variables
#
#_________________________________________________________________________________________
(TAX_ID, CHRM, CHR_START, CHR_STOP, STRAND, 
CONTIG, CTG_START, CTG_STOP, CTG_ORIENT, FEATURE_NAME, GENE_ID, 
FEATURE_TYPE, GROUP_LABEL, TRANSCRIPT, EVIDENCE_CODE) = range(15)
refseq_id_regex = re.compile("([^.]+(?:\..*)?)")# "NM_001011874.1"
def set_gene_and_transcript_variables (gene, values, locus, line_num, line):

    #
    # get cdna_refseq_id and transcript
    #   
    if values[FEATURE_TYPE] in set(("CDS", "RNA", "UTR")):
        matching_cdna_refseq_id = refseq_id_regex.search(values[TRANSCRIPT])
        if not matching_cdna_refseq_id:
            dump_and_die ("Transcript id not in the right format [%s]" % (values[TRANSCRIPT]), line_num, line)
        cdna_refseq_id = matching_cdna_refseq_id.group(1)

        # get transcript
        if not cdna_refseq_id in gene.transcripts:       # NM
            gene.transcripts[cdna_refseq_id] =t_transcript(cdna_refseq_id, gene) 
        transcript = gene.transcripts[cdna_refseq_id]



    #
    #   Save data for CDS, UTR, GENE, PSEUDO, RNA
    #
    if values[FEATURE_TYPE] == 'CDS':
        transcript.cds_loci.append(locus)        
        if values[FEATURE_NAME] != None  and values[FEATURE_NAME] != "" and values[FEATURE_NAME] != '-':

            matching_prot_refseq_id = refseq_id_regex.search(values[FEATURE_NAME])
            if not matching_prot_refseq_id:
                dump_and_die ("Refseq id not in the right format [%s]" % (values[FEATURE_NAME]), line_num, line)
            transcript.prot_refseq_id = matching_prot_refseq_id.group(1) 

    elif values[FEATURE_TYPE] == 'UTR':
        transcript.utr_loci.append(locus)


    elif values[FEATURE_TYPE] == 'GENE':
        gene.gene_name  = values[FEATURE_NAME]                   # Xkr4 or LOC100045372
        gene.evidence   = values[EVIDENCE_CODE]

    elif values[FEATURE_TYPE] == 'PSEUDO':
        if gene.gene_name != values[FEATURE_NAME]:
            status_log.die_error("pseudo gene %s has different gene name to gene %s on line #%d." %
                            (values[FEATURE_NAME], gene.gene_name, line_num))
        if gene.evidence != values[EVIDENCE_CODE]:
            status_log.die_error("pseudo gene %s has different evidence code to gene %s on line #%d." %
                            (values[EVIDENCE_CODE], gene.evidence, line_num))
        gene.pseudo_locus = locus

    elif values[FEATURE_TYPE] == 'RNA':
        transcript.evidence = values[EVIDENCE_CODE]
        transcript.rna_loci.append(locus)

#_________________________________________________________________________________________
# 
#   parse_exon_structure_file
#
#          parse seq_gene.md 
#
#_________________________________________________________________________________________
def parse_exon_structure_file (exonic_structures, gene_id_versions, exon_file, group_labels):
    """
    
        parses exon records from seq_gene.md into genes
        
    """
    # feature_name     gene_id            feature_type   transcript         evidence_code
    # ------------     -------            ------------   ----------         -------------
    #NP_001011874.1    GeneID:497097      CDS            NM_001011874.1     -
    #LOC100045372      GeneID:100045372   GENE           -                  protein;;
    #LOC100045372      GeneID:100045372   PSEUDO         -                  protein;;
    #XM_001472020.1    GeneID:100038975   RNA            XM_001472020.1     mRNA;identical;N
    #NM_001011874.1    GeneID:497097      UTR            NM_001011874.1     

    
    # progress dots begin
    if not options.quiet:
        sys.stderr.write("    [")
    cnt_genes_dots = [0,0, 1000]
    
    #
    #   regular expressions for parsing fields 
    #
    gene_id_regex = re.compile("GeneID:(\d+)")    # "GeneID:100131754"
        
    for line_num, line in enumerate( exon_file):
        
        #
        #   DEBUGG
        
        if options.debug_mode and line_num > 500000:
            break
        
        values = line.strip().split('\t')
        

        if values[GROUP_LABEL] not in group_labels:
            continue
            
    
        locus = t_locus(values[CHRM], int(values[CHR_START]) - 1, int(values[CHR_STOP]),
                        values[STRAND] == '+')
        
    
    
        # get gene_id
        matching_gene_id = gene_id_regex.search(values[GENE_ID])
        if not matching_gene_id:
            dump_and_die ("Gene id not in the right format [%s]" %(gene_id), line_num, line)
        gene_id = matching_gene_id.group(1)
        
        gene = create_or_retrieve_gene (values[FEATURE_TYPE],
                                        exonic_structures, gene_id_versions, 
                                        values[GROUP_LABEL], gene_id, locus, cnt_genes_dots)

                


        # use fields to update genes and transcripts
        set_gene_and_transcript_variables (gene, values, locus, line_num, line)
    
    
    # progress dots end
    if not options.quiet:
        sys.stderr.write("]\n")

    status_log.info("Finished parsing %d lines." % line_num)

    
    
    
    
    
    
#_________________________________________________________________________________________
# 
#   add_refseq_sequence
#          update transcripts with sequence and gis
#
#_________________________________________________________________________________________
from fasta_sequences_iterator import fasta_sequences_iterator
def add_refseq_sequence(seq_id_to_transcript, file_name, gi_field, seq_field, acc_field,
                        sequence_type_name):


    fasta_file_lines = open(file_name).readlines()
    unmatched_sequence_ids = set()

    if not options.quiet:
        sys.stderr.write("    [")
    cnt_genes_dots = [0,0, 1000]

    
            
    for seq_acc, sequence in fasta_sequences_iterator(fasta_file_lines):

        # get identifier and gi
        identifiers = seq_acc.split(' ')[0]
        identifiers = identifiers.split('|')
        versioned_seq_id = identifiers[3]
        #seq_id = versioned_seq_id.split('.')[0]
        seq_id = versioned_seq_id
        seq_gi = identifiers[1]
        
        # progress dots
        cnt_genes_dots[0] += 1
        if not options.quiet and (cnt_genes_dots[0]) % cnt_genes_dots[2] == 0:
            progress_dots(cnt_genes_dots)

        if not seq_id_to_transcript.has_key(seq_id):
            unmatched_sequence_ids.add(seq_id)
            continue

        # lookup transcript and set fields
        for transcript in seq_id_to_transcript[seq_id]:
            transcript.__dict__[gi_field] = seq_gi
            transcript.__dict__[seq_field] = sequence
            if acc_field:
                transcript.__dict__[acc_field] = seq_acc
        
    # end progress dots
    if not options.quiet:
        sys.stderr.write("]\n")

    status_log.info("%d %s sequences which did not have this group label" 
                    % (len(unmatched_sequence_ids), sequence_type_name))
        


#_________________________________________________________________________________________
# 
#   output_pseudos
#
#_________________________________________________________________________________________
def output_pseudos (all_group_label_genes):
    """
    write out list of pseudogenes
    """

    if not options.pseudos_with_cds_output_file:
        return
        
    #
    #   get pseudos
    # 
    pseudos_gene_ids_gene_types = list()
    for g in all_group_label_genes:
        g.get_pseudos_with_cds (pseudos_gene_ids_gene_types)
        
    #
    #   write out 
    #
    # make sure there is a directory ready for output
    prepare_directory_for_file(options.pseudos_with_cds_output_file)
    pseudos_with_cds_output = open(options.pseudos_with_cds_output_file, "w")
    for gene_id, gene_type in pseudos_gene_ids_gene_types:
        pseudos_with_cds_output.write( gene_type + "\t" + options.gene_accession_prefix + gene_id + "\n")
    pseudos_with_cds_output.close()
    
    
        


#_________________________________________________________________________________________
# 
#   output_coding_sequences_and_identifiers
#
#_________________________________________________________________________________________
def output_coding_sequences_and_identifiers (all_group_label_genes):
    """
        write out cdna and protein sequences for coding genes
        

    """

    # get accession and sequences
    acc_id_and_sequence = list()

    #
    #   get coding genes sequences
    # 
    for g in all_group_label_genes:
        g.get_acc_id_and_sequence (acc_id_and_sequence, 
                                    options.gene_accession_prefix, 
                                    options.prot_accession_prefix,
                                    options.cdna_accession_prefix,
                                    options.gis_for_cdna_prot_identifiers)

    #
    #   output cdna
    #
    if options.cdna_output_file:
        # make sure there is a directory ready for output
        prepare_directory_for_file(options.cdna_output_file)
        cdna_output = open(options.cdna_output_file, "w")
        for gene_id, acc, cdna_acc_id, prot_acc_id, cdna_seq, prot_seq in acc_id_and_sequence:
            cdna_output.write(">" + cdna_acc_id + " " + acc + "\n" + cdna_seq + "\n")
        cdna_output.close()

    #
    #   output prot
    #
    if options.prot_output_file:
        # make sure there is a directory ready for output
        prepare_directory_for_file(options.prot_output_file)
        prot_output = open(options.prot_output_file, "w")
        for gene_id, acc, cdna_acc_id, prot_acc_id, cdna_seq, prot_seq in acc_id_and_sequence:
            prot_output.write(">" + prot_acc_id + " " + acc + "\n" + prot_seq + "\n")
        prot_output.close()

    #
    #   output identifiers
    #
    if options.mapped_identifiers_output_file:
        # make sure there is a directory ready for output
        prepare_directory_for_file(options.mapped_identifiers_output_file)
        mapped_identifiers = open(options.mapped_identifiers_output_file, "w")
        for gene_id, acc, cdna_acc_id, prot_acc_id, cdna_seq, prot_seq in acc_id_and_sequence:
            mapped_identifiers.write("\t".join([gene_id, prot_acc_id, cdna_acc_id]) + "\n")
        mapped_identifiers.close()
        
        
def fields_to_tab_str (*fields):
    return "\t".join(map(str, fields))
#_________________________________________________________________________________________
# 
#   output_refseq
#
#_________________________________________________________________________________________
def output_refseq (all_group_label_genes):
    """
        print out all refseq for transcripts where the refseq and genomic sequence match 
        exactly
    """
    if not options.refseq_output_file:
        return
    #
    #   get coding genes refseq
    # 
    # make sure there is a directory ready for output
    prepare_directory_for_file(options.refseq_output_file)
    refseq_output = open(options.refseq_output_file, "w")
    gene_accession_prefix = options.gene_accession_prefix
    prot_accession_prefix = options.prot_accession_prefix
    cdna_accession_prefix = options.cdna_accession_prefix
    
    for g in all_group_label_genes:
        for t in g.transcripts.itervalues():
            if not len(t.cds_loci):
                continue
            genomic_match = t.refseq_genomic_match
            prot_refseq_id = t.prot_refseq_id
            if not prot_refseq_id:
                prot_refseq_id = ""
            cdna_refseq_id = t.cdna_refseq_id
            if not cdna_refseq_id:
                cdna_refseq_id = ""
            if not len(cdna_refseq_id) and not len(prot_refseq_id):
                continue
            if options.gis_for_cdna_prot_identifiers:
                prot_id = t.prot_gi
                cdna_id = t.cdna_gi
            else:
                prot_id = t.prot_refseq_id
                cdna_id = t.cdna_refseq_id
            refseq_output.write(fields_to_tab_str(
                                            gene_accession_prefix + g.gene_id, 
                                            prot_accession_prefix + prot_id,
                                            cdna_accession_prefix + cdna_id, 
                                            prot_refseq_id, cdna_refseq_id,
                                            genomic_match + 0) + "\n")

        





    #
    #   DEBUGG output in tabular format
    #
    #table_output = open("sequences.delete.this.table", "w")
    #for acc, cdna_acc_id, prot_acc_id, cdna_seq, prot_seq in acc_id_and_sequence:
    #    table_output.write("\t".join([acc, cdna_acc_id, prot_acc_id, cdna_seq, prot_seq]) + "\n")
    #
    #status_log.info("Written out sequence for %d CDS cdna and prot" % 
    #                        len(acc_id_and_sequence))


#_________________________________________________________________________________________
# 
#   output_coding_exonic_loci
#
#_________________________________________________________________________________________
def output_coding_exonic_loci (all_group_label_genes):
    """
        write out exonic_loci for all coding genes

    """
    if not options.exon_output_file:
        return
        
    # get accession and sequences
    exon_loci = list()

    #
    #   get coding genes
    # 
    for g in all_group_label_genes:
#        if g.gene_type != "coding":
#            continue
        g.get_coding_exon_list (exon_loci, options.gene_accession_prefix, 
                                            options.cdna_accession_prefix, 
                                            options.prot_accession_prefix,
                                            options.gis_for_cdna_prot_identifiers)

    # make sure there is a directory ready for output
    prepare_directory_for_file(options.exon_output_file)
    exon_output = open(options.exon_output_file, "w")
    for exon_data in exon_loci:
        exon_output.write("\t".join(exon_data) + "\n")
    exon_output.close()


#_________________________________________________________________________________________

#   dump_genes_to_files

#_________________________________________________________________________________________
def dump_genes_to_files(all_group_label_genes):

    def dump_exonic_structure(genes, output_file_name):
        # make sure there is a directory ready for output
        prepare_directory_for_file(output_file_name)
        cache_file = open(output_file_name + ".pic", 'wb')
        cPickle.dump(genes, cache_file, 1)
        cache_file.close()

        output_file = open(output_file_name, "w")
        for gene in genes:
            output_file.write(str(gene) + "\n")

    # 
    #   stops
    #
    stops_genes, mismatched_genes = get_genes_with_stops_mismatches (all_group_label_genes)
    if options.stops_genes_dump_file:
        status_log.info("Dumping genes with stops")        
        dump_exonic_structure(stops_genes,  options.stops_genes_dump_file)

    # 
    #   mismatches
    #
    if options.mismatched_genes_dump_file:
        status_log.info("Dumping genes with mismatches")        
        dump_exonic_structure(mismatched_genes,  options.mismatched_genes_dump_file)

    # 
    #   all genes
    #
    if options.all_genes_dump_file:
        status_log.info("Dumping all genes")        
        dump_exonic_structure(all_group_label_genes, options.all_genes_dump_file)
    # 
    #   gene_types
    #
    if options.all_gene_types_dump_file:
        status_log.info("Dumping gene types for all genes")        
        # make sure there is a directory ready for output
        prepare_directory_for_file(options.all_gene_types_dump_file)
        gene_type_dump = open(options.all_gene_types_dump_file, "w")
        for g in all_group_label_genes:
            gene_type_dump.write(options.gene_accession_prefix + g.gene_id + "\t" + g.gene_type + "\n")






    
#_________________________________________________________________________________________
#
#   print_transcript_count_with_duplicate_ids
#
#_________________________________________________________________________________________
def print_transcript_count_with_duplicate_ids (id_to_transcript, id_type_name):
    """
        takes dictionary of id to transcript
        counts transcripts which share the same identifier
    """
    cnt_duplicate_ids = defaultdict(lambda :0)
    for transcripts in id_to_transcript.itervalues():
        cnt_duplicate_ids[len(transcripts)] += 1

    duplicate_ids = [id for id, transcripts in id_to_transcript.iteritems() 
                            if len(transcripts) > 1]

    counts = sorted(cnt_duplicate_ids.keys())

    # 
    # no duplicates
    #
    if counts == [1]:
        return

    if duplicate_ids:
        status_log.info("Some transcripts shared %s identifiers" % id_type_name)
        for count in counts:
            #if count == 1:
            #    continue
            status_log.info("    %5d ids each found %5d times" %
                                (cnt_duplicate_ids[count], count))
        status_log.info("    Duplicate ids")
        for id in duplicate_ids:
            for transcript in id_to_transcript[id]:
                status_log.info("        %20s %20s" % (id, transcript.gene.gene_id))
        
    
    


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Main logic


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    


#
#    if exonic structure already exists in cache, use that
#        else parse seq_gene.md  
#
exonic_structures, gene_id_versions = get_exonic_structure_from_cache_file()



#
#   get list of all genes with requisite group label
#
all_group_label_genes =  get_group_labels_genes (exonic_structures, options.group_labels)





if count_PAR_gene_duplicates:
    status_log.info("%d gene_ids were shared between genes on PAR x/y chromosomes" % count_PAR_gene_duplicates)
    for gene_id in PAR_gene_duplicates:
        status_log.info("        " + options.gene_accession_prefix + gene_id)
if count_buggy_gene_duplicates:
    status_log.warning("%d gene_ids were duplicated probably due to a bug" % count_buggy_gene_duplicates)
    for gene_id in buggy_gene_duplicates:
        status_log.info("        " + options.gene_accession_prefix + gene_id)

        
#
#   add refseq sequence
# 
prot_refseq_id_to_transcript, cdna_refseq_id_to_transcript =   \
                        get_refseq_ids_to_transcript (all_group_label_genes)


print_transcript_count_with_duplicate_ids(prot_refseq_id_to_transcript, "prot refseq")
print_transcript_count_with_duplicate_ids(cdna_refseq_id_to_transcript, "cdna refseq")


status_log.info("Add refseq sequences...")
add_refseq_sequence(cdna_refseq_id_to_transcript,  
                    options.cdna_input_file, 
                    "cdna_gi", "cdna_refseq", "cdna_accession", "cdna")
add_refseq_sequence(prot_refseq_id_to_transcript,        
                    options.prot_input_file, 
                    "prot_gi", "prot_refseq", None, "peptide" )


#
#   check gis are unique
# 

prot_gi_to_transcript, cdna_gi_to_transcript =   \
                        get_gis_to_transcript (all_group_label_genes)

print_transcript_count_with_duplicate_ids(prot_gi_to_transcript, "prot gi")
print_transcript_count_with_duplicate_ids(cdna_gi_to_transcript, "cdna gi")
status_log.debug("Resolve duplicate GIs...")
resolve_duplicate_gis (prot_gi_to_transcript, "prot_gi")
resolve_duplicate_gis (cdna_gi_to_transcript, "cdna_gi")



#
#   add genomic sequence
# 

    
(exonic_coding_sequence, 
 coding_exonic_sequence_errors) = get_exonic_sequence_from_cache_file(all_group_label_genes)

   

#
#   add genomic sequence to genes
#  
add_genomic_coding_sequence (all_group_label_genes, exonic_coding_sequence, coding_exonic_sequence_errors)

status_log.debug("finalise_coding_sequence_data...")
for g  in all_group_label_genes:
    g.finalise_data()

#
#   output coding sequence statistics
#  
get_coding_sequence_statistics (all_group_label_genes)



#
#   output coding sequences
#  
output_coding_sequences_and_identifiers (all_group_label_genes)


#
#   output refseq associations
#  
output_refseq(all_group_label_genes)


#
#   output coding exonic loci
#  
output_coding_exonic_loci (all_group_label_genes)


#
#   output list of pseudo genes
#  
output_pseudos (all_group_label_genes)


#
#   dump genes for debug
#
dump_genes_to_files(all_group_label_genes)

# orderly shut down of status log
del status_log
del t_ncbi_gene_coding_structure.status_log


